#include "field.h" 

double sqr(double x)
{
	return(x*x); 
}
void field_tension(Field       *f_sigma_1_1,
				   Field       *f_sigma_1_2,
				   Field       *f_sigma_2_2, 
				   const Field *f_source1,
				   const Field *f_source2,
				   const Field *f_source3,
				   const Field *f_eta,
				   const Field *f_zeta,
				   const double K,
				   const double inv_dx,
				   const double inv_dy, 
				   const double inv_dx2, 
				   const double inv_dy2)
{
	CHECK_FIELD(f_sigma_1_1);
	CHECK_FIELD(f_sigma_1_2);
	CHECK_FIELD(f_sigma_2_2);
	CHECK_FIELD(f_source1);
	CHECK_FIELD(f_source2);
	CHECK_FIELD(f_source3);
	CHECK_FIELD(f_eta);
	CHECK_FIELD(f_zeta);
	CHECK_FIELDS(f_sigma_1_1, f_sigma_1_2);
	CHECK_FIELDS(f_sigma_1_1, f_sigma_2_2);
	CHECK_FIELDS(f_sigma_1_1, f_source1);
	CHECK_FIELDS(f_sigma_1_1, f_source2);
	CHECK_FIELDS(f_sigma_1_1, f_source3);
	CHECK_FIELDS(f_sigma_1_1, f_eta);
	CHECK_FIELDS(f_sigma_1_1, f_zeta);
	{
		const int nx = f_source1->nx;
		const int ny = f_source1->ny;
		const int nx1 = nx-1; 
		const int ny1 = ny-1; 
		double ** sigma11 = ARRAY_OF(f_sigma_1_1);
		double ** sigma12 = ARRAY_OF(f_sigma_1_2);
		double ** sigma22 = ARRAY_OF(f_sigma_2_2); 
		const double ** rho = CONST_ARRAY_OF(f_source1);
		const double ** rhu = CONST_ARRAY_OF(f_source2);
		const double ** rhv = CONST_ARRAY_OF(f_source3);
		const double ** eta1 = CONST_ARRAY_OF(f_eta);
		const double ** zeta1 = CONST_ARRAY_OF(f_zeta); 
		double fac_x = 0.5*inv_dx; 
		double fac_y = 0.5*inv_dy; 
		
		int i,j; 
		
		for(j=0;j<ny;++j)
        {
            const int jm = (j <= 0)   ? ny1 : j-1;
            const int jp = (j >= ny1) ? 0   : j+1;
            
            for(i=0;i<nx;++i)
            {
                const int im = (i <= 0)   ? nx1 : i-1;
                const int ip = (i >= nx1) ? 0   : i+1;
				
				sigma11[j][i] = 
				-0.5*K*(sqr(fac_x*(rho[j][ip]-rho[j][im]))-sqr(fac_y*(rho[jp][i]-rho[jm][i])))
				+K*rho[j][i]*(inv_dx2*(rho[j][ip]-2.0*rho[j][i]+rho[j][im])+inv_dy2*(rho[jp][i]-2.0*rho[j][i]+rho[jm][i]))+2.0*eta1[j][i]*inv_dx*(rhu[j][i]/(0.5*(rho[j][ip]+rho[j][i]))-rhu[j][im]/(0.5*(rho[j][i]+rho[j][im])))+(zeta1[j][i]-(2.0/3.0)*eta1[j][i])*(inv_dx*(rhu[j][i]/(0.5*(rho[j][ip]+rho[j][i]))-rhu[j][im]/(0.5*(rho[j][i]+rho[j][im])))+inv_dy*(rhv[j][i]/(0.5*(rho[jp][i]+rho[j][i]))-rhv[jm][i]/(0.5*(rho[j][i]+rho[jm][i]))));
				
				sigma12[j][i] = -K*(fac_x*(rho[j][ip]-rho[j][im])*fac_y*(rho[jp][i]-rho[jm][i]))+eta1[j][i]*(fac_x*(0.5*(rhv[j][ip]+rhv[jm][ip])/rho[j][ip]-0.5*(rhv[j][im]+rhv[jm][im])/rho[j][im])+fac_y*(0.5*(rhu[jp][i]+rhu[jp][im])/rho[jp][i]-0.5*(rhu[jm][i]+rhu[jm][im])/rho[jm][i]));
				
				sigma22[j][i] = 
				-0.5*K*(sqr((rho[jp][i]-rho[jm][i])*fac_y)-sqr((rho[j][ip]-rho[j][im])*fac_x))
				+K*rho[j][i]*((rho[j][ip]-2.0*rho[j][i]+rho[j][im])*inv_dx2+(rho[jp][i]-2.0*rho[j][i]+rho[jm][i])*inv_dy2)+2.0*eta1[j][i]*inv_dy*(rhv[j][i]/(0.5*(rho[jp][i]+rho[j][i]))-rhv[jm][i]/(0.5*(rho[j][i]+rho[jm][i])))+(zeta1[j][i]-(2.0/3.0)*eta1[j][i])*(inv_dx*(rhu[j][i]/(0.5*(rho[j][ip]+rho[j][i]))-rhu[j][im]/(0.5*(rho[j][i]+rho[j][im])))+inv_dy*(rhv[j][i]/(0.5*(rho[jp][i]+rho[j][i]))-rhv[jm][i]/(0.5*(rho[j][i]+rho[jm][i])))); 
			}
		}
	}
	
}

void field_tension_bis(Field       *f_sigma_1_1,
					   Field       *f_sigma_1_2,
					   Field       *f_sigma_2_2, 
					   const Field *f_source1,
					   const double K,
					   const double inv_dx,
					   const double inv_dy, 
					   const double inv_dx2, 
					   const double inv_dy2)
{
	CHECK_FIELD(f_sigma_1_1);
	CHECK_FIELD(f_sigma_1_2);
	CHECK_FIELD(f_sigma_2_2);
	CHECK_FIELD(f_source1);
	CHECK_FIELDS(f_sigma_1_1, f_sigma_1_2);
	CHECK_FIELDS(f_sigma_1_1, f_sigma_2_2);
	CHECK_FIELDS(f_sigma_1_1, f_source1);
	{
		const int nx = f_source1->nx;
		const int ny = f_source1->ny;
		const int nx1 = nx-1; 
		const int ny1 = ny-1; 
		double ** sigma11 = ARRAY_OF(f_sigma_1_1);
		double ** sigma12 = ARRAY_OF(f_sigma_1_2);
		double ** sigma22 = ARRAY_OF(f_sigma_2_2); 
		const double ** rho = CONST_ARRAY_OF(f_source1);
		
		double fac_x = 0.5*inv_dx; 
		double fac_y = 0.5*inv_dy; 
		
		int i,j; 
		
		for(j=0;j<ny;++j)
        {
            const int jm = (j <= 0)   ? ny1 : j-1;
            const int jp = (j >= ny1) ? 0   : j+1;
            
            for(i=0;i<nx;++i)
            {
                const int im = (i <= 0)   ? nx1 : i-1;
                const int ip = (i >= nx1) ? 0   : i+1;
				
				sigma11[j][i] = 
				-0.5*K*(sqr(fac_x*(rho[j][ip]-rho[j][im]))-sqr(fac_y*(rho[jp][i]-rho[jm][i])))
				+K*rho[j][i]*(inv_dx2*(rho[j][ip]-2.0*rho[j][i]+rho[j][im])+inv_dy2*(rho[jp][i]-2.0*rho[j][i]+rho[jm][i]));
				
				sigma12[j][i] = -K*(fac_x*(rho[j][ip]-rho[j][im])*fac_y*(rho[jp][i]-rho[jm][i]));
				
				sigma22[j][i] = 
				-0.5*K*(sqr((rho[jp][i]-rho[jm][i])*fac_y)-sqr((rho[j][ip]-rho[j][im])*fac_x))
				+K*rho[j][i]*((rho[j][ip]-2.0*rho[j][i]+rho[j][im])*inv_dx2+(rho[jp][i]-2.0*rho[j][i]+rho[jm][i])*inv_dy2);
				
				
			}
		}
	}
		
}